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Abstract 

Requirements are formulated for a reaction kinetics package to be useful for an as wide as possible circle of users and illustrated 

with examples using ReactionKinetics, a Mathematica based package. 

'Keywords: kinetics, mathematical modelling, dynamic simulation, computational chemistry, graphs of reactions, stochastic 

models 



1. Introduction 

■ In Part I of our paper we formulate the requirements for a 
reaction kinetics package to be useful for an as wide as possible 
circle of users. We try to answer the question: What should an 
ideal packag e know? 



In Part II ONagv et al 



201111 we enumerate the major prob- 



lems arising when writing and using such a package. 
' Throughout we try to illustrate everything with the present 
version of the package and the kinetic examples are mainly 
taken so as to join to the lectures presented at the workshop 
MaCKiE 2011 ( |http : //www . mack ie201 1 . uni-h d . de/| l. In 
hiany cases the examples will not show all the fine details nec- 
essary to be given when really using our program: a detailed 
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program documentation will be given later. Furthermore, triv- 
ial examples are chosen in some cases to transparently illustrate 
the functioning of the package. 

Also, there is no place to explain the theoretical background, 
we refer to the literature, and also to our future, more detailed 
work. 

2. History 

In the late sixties, early seventies of the last century, at the 
time when the computer became accessible for an ordinary sci- 
entist, kineticists immediately started to write codes for at least 
three different problems: for parsing large sets of reaction steps, 
solving induced kinetic differential equations and simulating 
the stochastic model o f chemical reactions, see an early review: 



GarfinkeletaL 



11197011 . 

Nowadays a few existing widely used program with multiple 
capabilities are e.g. 



CHEMKIN www . sandia . gov/chemkin/index . html 



XPPAUT www . math . pitt . edu/-bard/xpp/xpp . html 



KINALC garf ield. chem. elte .hu/Combustion/kinalc .htm 

July 18, 2011 



but there exists an almost infin ite number of t hem (w hich we are 



Tomlin et al. 



1 199711 ) with much 



going to review later; see also 
less capabilities. 

3. Coverage 



Our major governing principle is that the program should 
continuously include recent methods developed in most areas 
of reaction kinetics modeling. A conference like MaCKiE 201 1 
is an excellent occasion, a great help to learn newer methods. 

We show a few examples from the recent literature. 

3.1. Testing detailed balance 
3.1.1. Preliminaries and definition 



The long story starts possibly with 



Wegscheideillll901/2ll . but 



the first explicit formulation (without any fo rmula) of the prin 



Fowler and Milne 



ciple o f detailed balance has been given by 
1 192511 : in real thermodynamic equilibrium all the subprocesses 
(whatever they mean) should be in dynamic equilibrium sepa- 
rately in such a way that they do not stop but they proceed with 
the same velocity in both directions. 

One can say that this principle means that time is reversible at 
equilibrium, that is why the expression microscopic reversibil- 
ity is usually used as a synonym. 

The modern fo rmulation of the principle accepted by lUPAC 



Gold et al. 



1 199711 essentially means the same: "The principle of 
microscopic reversibility at equilibrium states that, in a system 
at equilibrium, any molecular process and the reverse of that 
process occur, on the average, at the same rate." 

Neither the above document nor the present authors assert 
that the principle should hold without any further assumptions; 
for us it is an important hypothesis the fulfilment of which 
should be checked individually in different models. 

Now let us consider a reaction consisting only of reversible 
steps: 

M 



m=\ 



M 



Y^a(m,p)X(m)f^Yj/3(m,p)Xim), (p = 1,2, . . . ,P) (1) 



where we have P pairs of reaction steps, a, ft are the matrices 
with molecularities as elements, their diff'erence is the stoichio- 
metric matrix, and kp (p = -P, . . . , -2, - 1 , 1 , 2, . . . , f ) are the 
reaction rate coefficients. Now, the induced kinetic differential 
equation of the above reaction is as follows: 

c(0 = 
p 

J]m; P) - «(% P)) [kp^tf^'''' - k-pcitf^'''^^ (2) 

p=\ 

where c(0 e M.'^ describes the concentrations of the species at 

time t. 

Definition 1 Assume that c* e (R+)*' and it makes the right- 
hand side of (|2|l to be equal to zero, furthermore the condition 



kp(er-^"^^k^„(ef^'"^ 



(3) 



is fulfilled for every p - 1,2 . . .,P. Then we say that reaction 
([T]i is detailed balanced at the stationary point c*. If the reac- 
tion is detailed balanced at all of its positive stationary points, 
then it is detailed balanced. 

3.1.2. Methods to test detailed balance 

The most naive way to check conditions ^ is the direct one: 
having found one (or all of) the positive stationary points sub- 
stitute it (them) into ^ and see if it holds or not. Or, it may 
be enough to solve equations (O for positive c*'s to get all the 
candidates in which reaction ([U is detailed balanced. However 
a nonlinear system of equations is required to be solved, which 
makes the problem difficult. 

Fortunate ly, there is a mor e elegant way to look at the prob- 



lem, namelv lFeinberd 0198911 has proved that the following two 
conditions are necessary and sufficient for the reaction ([Ij be 
detailed balanced. 

Condition 1 (circuit conditions) Suppose that we have cho- 
sen an arbitrary spanning forest for the Feinberg-Horn-Jackson 
graph of reaction ([TJ. It is possible to find a set of P - A^ -H L 
independent circuits, where A^ denotes the number of vertices 
{complexes): 



N 



{Qf(.,r);r= 1,2, ...,/?) U {)8(., r); r = 1,2, 



,^!l, 



whereas L is the number of connected components (linkage 
classes) of the Feinberg-Horn-Jackson graph (in which all the 
complexes are written down exactly once and they are con- 
nected with reaction arrows). For each of these circuits we write 
an equation which asserts that the product of the reaction rate 
coefficients in the clockwise direction and counterclockwise di- 
rection is equal. Thus we have P - N + L equations for the 
reaction rate coefficients. 

Before formulating Condition 2 we need another important def- 
inition. 

Definition 2 The deficiency is the number of complexes A^ mi- 
nus the number L of connected components of the Feinberg- 
Horn-Jackson graph and the number S of independent reaction 
steps (or, the rank of the stoichiometric matrix). 

Condition 2 (spanning forest conditions) Assume that reaction 
([Ij shows deficiency 6. Furthermore, assume that the edges of 
an arbitrarily selected the spanning forest F has been given an 
orientation. Then there exists 6 independent non-trivial solu- 
tions to the vector equation 

2 (fi{;j)-a{:,i))a{i,j)^Q 

for the a{i, j) numbers. With these solutions one can construct 
the spanning forest conditions which are 

(.i,i)eF (i,j)eF 

where k(i, y)'s are the corresponding reaction rate coefficients 
(associated to the edge (/, j)). 

3.1.3. Applications 

Here we present only one example. Take the model of 
Wegscheider, i.e. GetProblem["Wegscheider"] to obtain 

{"A" ^ "B", 2 "A" ^ "A" + "B"). 

Now 

DetailedBalanced ["Wegscheider" ,ki ,k_i ,k2,k_2] 
results in 



DetailedBalanced: :nocycle: The FHJ graph of the given formal 
mechanism has no cycle, so we are given only the spanning forest 
condition(s). 

{k_i*k2 == k_2*ki} 

providing the condition for the reaction to be detailed balanced. 
There are several applications for e.g. in chirality, models of 
ion channels, combustion theory in progress. 

3.2. Absolute concentration robustness 

One might look for sufficient conditions to ensure the inde- 
pendence of a stationary concentration coordinate from outer 
conditions: does it only depe nd on reaction rate coefficien ts and 



on no initial concentrations? IShinar and Feinberg , 



201C I. 



iA + B-^2B,B^AL 



3.2.1. Motivation 

A simple example follows. Let us consider the reaction 

rl := 1^ 

where A may be interpreted as the inactive form of a protein 
and B as its active form, e. g. trypsinogen and trypsin re- 
spectively. (Let us remark in passing that this model is an ir- 
reversible version of the Wegscheider model above.) The in- 
duced kinetic differential equation of this reaction is given by 
DeterministicModel [rl, {ki , k2)]] as 

a'(t) - -kia(t)bit) + k2b(t) b'(t) - kia(t)b(t) - k2b(t). 

Calculation of the stationary points (under the condition that 
a(0) = flo, ^(0) = bo) proceeds quite naturally: 
StationaryPoints [rl , {ki , k2), {ao, bo}, 
Conditions -> jao + bo > k2/ki), 
Positivity -> True] 

The stationary concentration fl*of A, being j^, is always pos- 
itive (together with b* - ao + bo - j^, at least if ao + bo > p), 
and does not depend on total mass: the reaction shows absolute 
concentration robustness with respect to A. 

However, in the reaction A ^ B neither A, nor B has this 
property, both coordinates of the stationary concentration vec- 
tor do depend on oq + bo (as one can easily find it without any 
program): «*= A:: 1^, fe'=^-i|^. 



[±r 



-kr 



more complicated cases. Nevertheless one should be cautious. 



Figure 1: The Feinberg-Horn-Jackson graph of the reaction r\ ob- 
tained by ShowFHJGraph[rl,Style[#,Red,14]&/(a|ki , k2|, 
VertexLabeling->True,DirectedEdges->True] 

The question is what kind of general conditions can be given 
to assure absolute concentration robustne ss. An easy /to/check 



set of conditions has been formulated by 



I 201C I 



Shinar and Feinberg 



3.2.2. Conditions 

Theorem 1 Suppose that a reaction endowed with mass action 
kinetics has a positive stationary point, and its deficiency is one. 
Then, if there exist two complexes lying in nonterminal strong 
components of the Feinberg-Horn-Jackson graph which only 
differ in a single species, then the reaction shows absolute con- 
centration robustness with respect to this species. 

The notions in the theorem are widely used and can be 
found in the paper cited above. In the case of the simple 
example above: ReactionsData[rl] ["deficiency"] gives 
5 = 4-2-1 = 1. 

We have shown earlier that there exists a positive stationary 
point. 

To see that there exist two nonterminal complexes only dif- 
fering in a single species we draw the Feinberg-Horn-Jackson 
graph of the reaction, and gladly observe that (A -H B) - A = B. 
(To make it clear, complex A H- 2B and B differs in both A and 
B.) 



The last condition is also checked by our program, it would be 
really hard to test this property in a large reaction set by hand. 

Here is how to learn which species of the given reaction are 
absolutely robust: AbsoluteRobustness [rl] gives (A). 

Let us remark that the program is also capable to calculate the 
stationary concentrations both numerically and symbolically in 



Remark 1 The theoretical and computational problem lies in 
showing that a positive stationary point exists. Why? Be- 
cause if the initial value problem to describe a reaction is 
c' = foe, c(0) = Co then the chemically meaningful 
stationary point c* 

1. should obey f(c*) = 0, 

2. should be nonnegative, 

3. should be on the level sets of the linear first integrals, and 

4. should be on the level sets of the nonlinear first integrals. 

The situation is that one cannot in general find (global) nonlin- 
ear first integrals. 

Let us consider reaction rl in detail. The solutions to f (c*) = 



in this case are 







(a* e 



and 



b* 



(b* 6 M). 



These solutions are nonnegative if a* > and b* > 0, respec- 
tively. They are on the level curves of the linear first integral 



i//(p, q) :- p + q if and only if < ciq + bo 



No nonlinear 



global first integral — independent from i/r — exists. 

Let us return to detailed balancing. Certainly, the con- 
ditions were not found to treat the trivial examples above. 
Their role will be clearer if we cite another — this time more 
complicated — example from the paper. 

Let us consider a model where ATP is the cofactor in the 
osmoregulation system of Escherichia coli: EnvZ-OmpR. 

r5 = {X ^ XT ^ Xp,Xp +Y <^ X,,Y ^ X -i- Y^, 
XT + Yp^ XTYp ^XT + Y); 

The existence of a positive stationary point can be proved either 
numerically, or symbolically using StationaryPoints [r5, 
Positivity -> True] and — waiting for a very long time. 

The result given by React ionsData[r5] ["deficiency"] 
is5 = n-l-s = 9-3-5 = l. 

To see that there exist two nonterminal complexes only dif- 
fering in a single species we draw the Feinberg-Horn-Jackson 
graph of the reaction (Figure [3.2.2b and find that (XT + Y p) — 



E- 



\i'^ 



Subsequent maxima 



Kp + X I 



Figure 2: The Feinberg-Horn-Jackson graph of reaction r5 obtained by 
ShowFHJGraph [r5 , VertexLabeliiig->True , DlrectedEdges->True] 
with proper complex colourings (purple ones are the terminal complexes). 

XT = Yp, therefore the reaction is absolutely robust with re- 
spect to Yp, the phosphorylated form of the response-regulator, 
a species of crucial importance in the reaction. 



3.2.3. The conditions are only sufficient but not necessary 

The fact that the conditions are only sufficient but not 
necessary can also be s een on the examples given by 



I Shinar and Feinberg , 



2010(1 . see also the Supplementary mate- 



rial. Here we give another ex ample, the o ne shown by Professor 



Ross in his lecture (see also JRoss , 



20081 p. 21.3611: 



;'" = 
^ Xi ^ X2 ^ X3 ^ X4 ^ Xs ^ Xfi ^ Xv ^ Xg ^ 

(4) 

The result given by ReactionsData[jr] ["deficiency"] 
is (5 = n - I - s - 9-1-8 = 0, (no wonder, 
jr is a compartmental system), the theorem cannot be applied, 
although all the stationary concentrations are only dependent 
on the reaction rate coefficients. Expression 

StationaryPoint [jr, 

(ko,ki,k_i,k2,k_2,k3,k-3,k4,k_4, 

ks , k_5 , ke , k_6 , k? , k_7 , kg} , 

{ci [0] , C2 [0] , C3 [0] , C4 [0] , C5 [0] , C6 [0] , c^ [0] , cg [0] }] 

gives the result symbolically, what we do not reproduce here 
verbatim, because of the length of the formulae. Here is an 




Figure 3: Return to the stationary state after a perturbation of the first species 
in the reaction jj) with the reaction rate coefficients J5)- 



(/= 1,2, 



,8) 



equivalent symbolic form: 

the proof and generalization of which is left to the reader. 

Once here, we can also qualitatively reproduce Fig. 3 of the 
paper bv lRossI 1200811 . Let us use the following set of reaction 



rate coefficients: 

rrc = {0.1,2,0.1,8,5,3,0.4,1,1,6,0.5,4,2,10,1,1} 

(5) 
and let the initial concentrations be the same as the stationary 
concentrations except that we add 100 to the initial value of the 
initial concentration of Xi. Than the return of the perturbed 
concentrations — which was calculated by 

Concentrations [jr, rrc, ini+lOOUnitVector [8, 1] , 
(0,4}]] 

— to the (asymptotically stable) original stationary state is as 
seen on our Fig. 13.2.31 



More examples and many further interesting details can be 
found in the mentioned Science paper. 

3.3. Improved methods of stochastic simulation 
As in Part II of our paper 



Nagv et al. 



1 201111 we have exposed 



we have built in all the relevant direct and approximate methods 



number of species 

i 



Species: X, Y 




Figure 4: Stochastic Lotka-Volteira model with reaction rate coefficients ki = 
\, k2 = 1/1000, ^3 = 1 and initial conditions xq = 600, yo = 400, using tau- 
leaping methods and compiled functions, where 2.157 CPU time units were 
elapsed. 

as well as explicit and implicit ones into our program package. 
For a long while several fancy improvements are known and are 
involved, here we intend to mention only two ways. 

The first one is more theoretical and relies on numerical 
techniques which have already been applied successfully when 
solving ordinary differentia l equations. We mentioned a few of 



them in 



Nagv et al 



11201111 but further ideas can and will im- 
prove the approximation methods of stochastic simulation of 
chemical reactions. 

The other way is based on recent developments and pro- 
gramming tricks. The latest versions of Mathematica enable 
us to use compiled functions in a more efficient and delicate 
way which are e.g. thousands of times more effective in func- 
tion evaluation compared to the "usual" evaluations (see the 
function Compile). Another direction of recent developments 
concerns parallel computing (see the functions Parallelize, 
ParallelMap, etc.), and the use of GPUs (see the functions 
CUDALink and OpenCLLink, etc.). 

Here stands an example: the Lotka-Volterra model is 



k] k^ k^ 

X -U 2X, X H- Y — * 2Y, Y ^ 0, 



and see Figure 1331 for the simulation results. 



3.4. Methods of reaction generation and decomposition 

A fundamental problem of stoichiometry is the decomposi- 
tion of overall reactions into elementary steps. The definition 
of "elementary" (or "simple") reactions varies; we call here a 
step elementary if it is of order at most two. Ideally, a reac- 
tion kinetics package should not only provide the user with the 
list of possible decompositions (or at least a large number of 
decompositions, as their number may be huge, or even infi- 
nite), but it should also assist the user in the identification of 
elementary steps themselves, or even in the generation of pos- 
sible intermediate species that may take part in an elementary 
step. Clearly, the complete automation of this three-step pro- 
cess is a formidable task, as a large amount of domain specific 
expertise and data must be incorporated into each step. A more 
reasonable goal is to provide the user with a generous list of so- 
lutions, all of which adhere to the basic conservation laws (e.g., 
all elementary steps must conserve the number of atoms of each 
element and total charge) and other combinatorial constraints 
imposed by the problem; and leave all further processing to the 
user. At this abstract level the generation of elementary steps 
with given reactants and the generation of decompositions of a 
given overall reaction are essentially equivalent, as we shall see 
immediately. 

Suppose that the species are made of D different elements, 
and assign a (D -i- l)-dimensional vector a,„ to species X{m) 
where the first D components are the quantities of each con- 
stituent, and the last component is electric charge. A similar 
vector can be assigned analogously to every complex as well. 
Now, a reaction 

M M 

2_^amX{m) — > 2_jP»'^^^) 

m= 1 m= 1 

obeys the laws of atomic and charge balance if and only if the 
vectors a^ describing the atomic structure of the species satisfy 
the linear system of equations 

M M 

m= 1 m=l 

All combinatorially feasible elementary reactions can be gener- 
ated by expressing each possible reactant complex (correspond- 



ing to a vector 'Zm^i^m^m satisfying 1 < Y,m<^m < 2) as a 
linear combination of the vectors a^, with nonnegative integer 
coefficients. Hence, we are looking for the nonnegative integer 
solutions X of a system of linear equations 

Ax = b, 

where A is the atomic matrix, with columns ai through slm (M 
denotes the number of species), and b is a vector corresponding 
to the reactants. With M species the elementary reactions are 
obtained by solving 2M + (^| such systems. 

Similarly, if M is the number of species, an M-dimensional 
vector is associated to every reaction involving these species. 
The mth component of the vector shows the change in the quan- 
tity of species X{m) in the reaction; if a species is only a reac- 
tant, its component is negative, if it is only a product, the co- 
efficient is positive, if it occurs on both sides, the coefficient 
might be either positive or negative, or even zero. The matrix / 
with column vectors corresponding to the elementary reactions 
is the stoichiometric matrix of the mechanism. A combination 
of elementary steps with coefficients x is a decomposition of 
the overall reaction if and only if 

■yx-w, 

where w is the vector associated with the overall reaction. 
Again, we are primarily interested in nonnegative integer so- 
lutions X, though nonnegative rational solutions also can be in- 
terpreted as decompositions of the overall reaction. 

While the complexity of this problem is very high (indeed, 
deciding whether there exists even a single solution is NP-hard), 
research in these fields have yielded several relatively practi- 
cal algorithms and heuristics for its solution. Further cross- 
fertiUzation between fields in this problem is essential, but one 
must also keep in mind that different applications produce prob- 
lem instances with strikingly different characteristics, which in 
tur n call for different m etho ds. Indeed, in our inve stigations 



in BKovacset al. . 



200411 and BPappand Vizvari , 



200611 we have 



found that entirely different algorithms are the most effective in 
elementary step generation and in the generation of the decom- 
positions. In particular, note the following obvious differences: 



• in elementary step generation multiple relatively small 
systems need to be solved, all of which share the coef- 
ficient matrix, whereas in decomposition a single large- 
scale system is solved; 

• in elementary step generation all but one row consist of 
nonnegative numbers only, in decomposition every row 
and column may contain entries of different sign; 

• in elementary step generation the number of solutions is 
always finite, in decomposition it is often infinite. 

Elaborating on the the last point, if a sequence of steps forms 
a cycle in which no species are generated or consumed, it can 
be added to any decomposition an arbitrary number of times. 
A stra ightforward application of Dickson's lemma 
1 191311 show is that the converse is also true: 



Dickson 



Lemma 1 The number of decompositions is finite if and only if 
the elementary reactions cannot form a cycle. Furthermore, the 
number of decompositions that do not contain a cycle is always 
finite. 

The lemma motivates two further problems related to decom- 
positions: first, one needs to be able to decide whether the el- 
ementary reactions can form cycles or not, which amounts to 
the solution of a linear programming problem. Second, if cy- 
cles do exist, the goal changes from generating all decompo- 
sitions to generating all minimal (that is, cycle-free) decompo- 
sitions, and all minimal cycles (that is, cycles that cannot be 
expressed as a sum of two cycles). The latter is equivalent to 
the problem of generat i ng all minimal P-invariants of a Petri 



net JMartinez and Silva , 



1982ll . 



Another approach to handle the explosion in the number of 
solutions is to restrict the attention to the simplest decomposi- 
tions and cycles, those that consist of a small number of steps. 

Returning to our discussion of basic requirements for a gen- 
eral purpose reaction kinetics package, it has become clear to us 
that our package must implement a number of methods for the 
problems of reaction generation and decomposition, along with 



Ag+ 


Ag2- 


H+ 


SO4 


S02- 


s2or 


C2O2- 


Ag(C204) 


OH 


H20 


co- 


02 


HO2 


H2O2 


O2CO2 


C02 



Table 1 : Species of the oxalate-persulphate-silver oscillator 

an algorithm selection heuristic to choose the one most suited 
for the problem at hand. 

Currently three algorithms are implemented for the solution 
of these problems, including one of our own developed with 
very large-scale problems in mind. Additionally, we provide 
two preprocessing methods to identify the steps that must take 
part in every decomposition, and the ones that cannot take part 
in any. (This greatly reduces the computation time for some of 
the algorithms.) We also provide a heuristic that generates a 
typically large number of decompositions much faster than the 
rigorous algorithms, but without any guarantee that it finds all 
minimal decompositions. 



3.4.1. Example 

We conclude this section with an example for elemen- 
tary step gene ration. The oxalate-persulphate-silver oscillator 



I Clarke . 



199211 involves 16 species, shown in Table [T] these can 



form 2- 16 + (2) = 152 complexes that may be the reactants 
of an elementary step. The species consist of five atomic con- 
stituents, and have charge, hence they can be represented by a 
6x16 atomic matrix, generated by the AtomMatrix command 
of our package. From this matrix the Elementary-Reactions 
command generates all elementary steps involving them. In this 
tiny example each algorithm currently implemented proved to 
be rather efficient; the number of combinatorially possible ele- 
mentary steps is 89. 

In Part II of our paper we shall revisit this problem, and 
provide a detailed example of obtaining a decomposition for 
a complex overall reaction. 



4. Applications 

Although there are big differences in the different fields of 
chemical kinetics, e.g. the role of thermodynamic data is more 
important in combustion, much less important in inorganic 
chemistry, still we hope the package will be used for more and 
more realistic applications in many fields. Here we only men- 
tion a few of our previous applications with the ancestors of the 
present package. 



Reactions on a surface, reactions in plasma: 



1 1974 1. 



Sipos et al 



Enzyme kinetics: 
Signal transduction: 
Ion channels: 



TothI 112002 1. 



Toth and RosparsI ll2005ll . 



Nagv et al. 



1 2009 1. 



5. Formats 

A useful program should be compatible with usual formats 
such as CHEMKIN, Prime, SBML etc. What we can do at 
the moment is that we can read in data of different formats 
but not in an automatic way. The main method to transform 
data in different forms is to use pattern matching, including 
if necessary or useful such tools as RegularExpression. 
Suppose we have the file hydrox.dat from the website 
http : //www . math . bme . hu/ -nagyall describing a model of 
hydrogen oxidation. Then, we can import and transform it like 
this: 

Import ["hydrox.dat"] /. {x , "SPECIES", 

Shortest [y ], "END", z } :> y The problem is 

obviously similar to, or may be considered to be part of 
parsing. 

6. Speed and Accuracy 

It is important to be able to handle a reaction during a rea- 
sonable time interval even if it consists of a large number of 
reaction steps and species. Similarly, the accuracy should also 
be enough for comparisons with measurements. 



Let us mention a few more specific problems where the num- 
ber of species and of reaction steps can be really high — beyond 
the well known areas of atmospheric chemistry, metabolism and 
combustion. One might wish to treat reaction chromatogra- 
phy in such a way that one divides a column into thousands 
of plates, assumes the same reaction on each of the plates and 
also assumes (linear) diffusion between the plates. This is the 
problem mentioned by Prof. Trapp in his lecture at the con- 
ference. Let us also ment ion here the less known paper by 



publications) within the same framework. In num erical com- 



Shapiro and Horn 



Ill979albll which gives a qualitative treatment 
of such systems using the tools of Chemical Reaction Network 
Theory. 

Or, one might describe the transformation taking place 
among thousands of polymers with different molecular 
weight, as mentione d in the poster by P. van Steenberge 
van Steenberge et al. . 

Another problem might be the treatment of molecules sitting 
at different energy levels of which one might have several mil- 
Uons. A possible first step to treat such a system might be to 
measure the time needed to solve the induced kinetic differen- 
tial equation of a model as a function of size as follows. 

lendvayCnJ := Table [X,- <-^ X,+ i, {i, 1, n - 1)] 
Concentrations [lendvay [1000] , Array[N[#] &, 
1998], Array [0.1 N [#] &, 1000], {0, 1)] // 
Timing 

It turned out that such a model can be solved in 30 
seconds using Mathematica without our package, without 
Parallelize and without using GPU and similar tools, i. e. 
one can say, only applying Low Performance Computing. This 
is quite a promising start in this direction. 

7. Language: Mathematica, what else? 

We do not want to enter into an infinitely long discussion 
about the advantages and disadvantages of mathematical pro- 
gram packages, we are only making a few remarks, which are 
possibly acceptable by the majority of our readers. 

Mathematica is capable of symbolic and numeric calcula- 
tions, creating graphics (in fact, even whole presentations and 



puta tions it is not worse than any other program llWeissteinI : 



Inc. 11 . It uses as many cores/processors as you have in paral- 
lel, it implements OpenCL and CUDA to use GPUs, you can 
ask it to calculate the C form of a compiled function if you 
need, etc.) Using a symbolic-numeric mathematical program- 
ming language such as Mathematica is extremely helpful in cre- 
ating very transparent programs; code that can even be read as 
"pseudo-code" for those readers who are most reluctant to get 
acquainted with the software. 

In a meticu l ous a nalysis of the consecutive reaction 



Yablonskv et al. 



1 201011 have found a symbolically expressed 
necessary and sufficient condition for the three concentration 
time curves to h ave a single point of intersection (see also 



I Toth and SimonL page 341]). As a final application of our 
package we sow how to solve this (and also similar, symbol- 
ically untractable problems) numerically. 



Intersecting concentration time curves 




/.-,- 
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Figure 5: Concentration time curves of tlie consecutive reaction 
obtained by Manipulate [ Plot [Evaluate [ ReplaceAll 00 

Concentrations [cons, k, 1, 0, 0, 0, 3]], t, 0, 3] 



8. Outlook 

Having collected so much requirements with illustrations in 
the second part we are going to turn to the problems and diffi- 



culties when writing and using a program package pretending 
to solve so many problems of reaction kinetics. 
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